%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Ramsey fringes with Bloch equations including relaxation.
%  n  = (n11 n22 n12 n21) .
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Bloch
    figure(1);					clf;
    col='rgbmcky';
    
    O  =            1;                                                      % Ramsey pulse power
    T  =       pi/2/O;                                                      % Ramsey pulse duration
    tau=            8;                                                      % time between pulses
    G  =         .1*O;                                                      % decay rate
    D  =[-1:.002:1]*8;                                                      % Ramsey pulse detuning
    psi0=[1 0]';
    n0 =[1 0 0 0]';

% evolution on the Bloch sphere

    dt =.2;
    kMx=15;
    for k=1:kMx
        n =n0;

        for tt=0:dt:T+tau+T
            d = 0;                          o = O;
            if tt<T+tau
                d =-.04*(k-ceil(kMx/2));   o = 0;
            end
            if tt<T
                d = 0;                      o = O;
            end

            ns(1)=sqrt(2)*imag(n(3));
            ns(2)=sqrt(2)*real(n(3));
            ns(3)=real(n(2)-n(1));
            n = Bloch2(n,dt,0,o,d);
            nt(1)=sqrt(2)*imag(n(3));
            nt(2)=sqrt(2)*real(n(3));
            nt(3)=real(n(2)-n(1));

            subplot(231),plot3([ns(1) nt(1)],[ns(2) nt(2)],[ns(3) nt(3)],'b','MarkerSize',10,'LineWidth',1.2)
            set(gca,'xLim',[-1 1],'yLim',[-1 1],'zLim',[-1 1],'FontSize',18,'LineWidth',1.2)
            xlabel('\fontsize{18}Im {\it\rho}_{12}')
            ylabel('\fontsize{18}Re {\it\rho}_{12}')
            zlabel('\fontsize{18}2{\it\rho}_{22}-1')
            text(1,1.5,1.3,'\fontsize{18}(a)')
            view([-110 20])
            hold on
            pause(.01);
        end
    end
    hold off
    grid

% Schrödinger equation

    for dd=1:length(D)
        psi=Schrd2(psi0,T,O,D(dd));
        psi=Schrd2(psi,tau,0,D(dd));
        psi=Schrd2(psi,T,O,D(dd));
        nt(dd)=abs(psi(2)).^2;
    end
    subplot(232),plot(D,nt,'LineWidth',.8)
    set(gca,'xLim',[-1 1]*max(D),'FontSize',18,'LineWidth',1.5)
    xlabel('\fontsize{18}\Delta_{12}')
    ylabel('\fontsize{18}{\it\rho}_{22}')
    text(-12,1,'\fontsize{18}(b)')
    hold on

% Bloch equation

    for dd=1:length(D)
        n=Bloch2(n0,T,G,O,D(dd));
        n=Bloch2(n,tau,G,0,D(dd));
        n=Bloch2(n,T,G,O,D(dd));
        nt(dd)=n(2);
    end
    subplot(232),plot(D,nt,'r','LineWidth',1.5)
    set(gca,'FontSize',18,'LineWidth',1.5)
    xlabel('\fontsize{18}\Delta_{12}')
    ylabel('\fontsize{18}{\it\rho}_{22}')
    pause(.01);

% Bloch equation 3D
    
    Tau=[0:.05:5]*tau;
    for tt=1:length(Tau)
        for dd=1:length(D)
            n=Bloch2(n0,T,0,O,.3*D(dd));
            n=Bloch2(n,Tau(tt),0,0,.3*D(dd));
            n=Bloch2(n,T,0,O,.3*D(dd));
            nt(tt,dd)=real(n(2));
        end
    end
    subplot(233),mesh(D,Tau,nt)
    set(gca,'xLim',[-1 1]*max(D),'yLim',[0 40],'FontSize',18,'LineWidth',1.5)
    xlabel('\fontsize{18}\Delta_{12}')
    ylabel('\fontsize{18}{\itT}')
    zlabel('\fontsize{18}{\it\rho}_{22}')
    view([-20 60])
    text(-12,40,1.5,'\fontsize{18}(c)')
    pause(.01);
end

function psi = Schrd2(psi0,t,o,d)
    H =[0    o/2 
        o/2   d ];
    psi=expm(1i*H*t)*psi0;
end

function n = Bloch2(n0,t,g,o,d)
    I =     1i/2;

    M=[ 0      g     I*o   -I*o 
        0     -g	-I*o    I*o 
        I*o  -I*o -g/2-1i*d  0   
       -I*o   I*o     0   -g/2+1i*d];
    n=expm(M*t)*n0;
end

